      subroutine preclinsol(x,b,Nt,PP, Q, R, S, T,
     >                      UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb,
     >                      Ycomm,va)
C
C     June 8 1998, Auke van der Ploeg, MARIN.
C
C     This subroutine solves the system of 
C     linear equations Ax=b, in which A is the Jacobian matrix that
C     corresponds to L1b xi-constant planes. 
C     The routine constructs a coupled incomplete LU-decomposition 
C     which is used as a preconditioner for GMRES.
C
C
C  Parnassos uses 4*4-block structure of the matrices.
C
C   row 1. momentum equation in main stream direction (ksi i-coordinate)
C   row 2. continuity equation
C   row 3. momentum equation in girthwise direction (zeta k-coordinate)
C   row 4. momentum equation in wall-normal direction (eta j-coordinate)
C
C   First column of 4*4 block multiplies U1 (velocity in x-direction).
C   Second column of 4*4 block multiplies U2 (velocity in y-direction).
C   Third column of 4*4 block multiplies U3 (velocity in z-direction).
C   Fourth column of 4*4 block multiplies P  (pressure).
C
C  Matrix entry locations
C
C               T          k
C               |          |
C          PP - Q - R      +- j
C               |   
C               S 
C
C  DOWNSTRa is coefficient in downstream direction (i+1)
C  DOWNSTRb is coefficient in downstream direction (i+2)
C  UPSTRa is coefficient in upstream direction (i-1)
C  UPSTRb is coefficient in upstream direction (i-2)
C=======================================================================
C         Description of Arguments
C=======================================================================
C 
C     x:  on input, current approximation to the solution.
C     x:  on output, the correction for the solution is stored in x.
C     b:      i/o right-hand side, on output overwritten by 
C                 preconditioned right-hand side.
C     Nt: size of A.
C     tolgmr: GMRES is stopped when
C             the 2-norm or the absolute norm of the residual has
C             decreased by a factor tolgmr 
C
      include 'dataparnas.incl'
      include 'dataparal.incl'
      include 'datamat.incl'
      include 'datava.incl'
      external matvec,ope
      INTEGER Nt,cnt,nb1
      REAL*8 x(Nt),b(Nt),dnrm2own,tol2
C
C     e1,e2,e4 and e5: coefficients required for the preconditioner:
C
      REAL*8,allocatable,DIMENSION(:,:,:,:,:,:):: e1,e5,e2,e4
C
C=======================================================================
C        PARAMETERS for GMRES
C=======================================================================
C
      INTEGER restrtgmr, maxliniter, iabs, ioutgmr,lwork,lmax
      LOGICAL nonzerox,rightprecon,knoll
      PARAMETER (lmax=2, knoll=.FALSE., rightprecon=.FALSE.)
      REAL*8,allocatable,dimension(:)::  work
      REAL*8 tol1,R2normprecloc,R2normloc,Fraction,ff
      INTEGER iflag, liniter, restrt

      restrtgmr=maxlingmr
      maxliniter=2*maxlingmr
      lwork=restrtgmr+1

!$OMP BARRIER
!$OMP FLUSH
!
!  Allocate local arrays for the preconditioner 
!
      ALLOCATE(work(lwork*Nt),
     >         e1(N,N,NY,NZ,1:1,L1b),e4(N,N,NY,NZ,0:0,L1b),
     >         e5(N,N,NY,NZ,1:1,L1b),e2(N,N,NY,NZ,0:0,L1b),stat=err)
!$OMP BARRIER
!$OMP FLUSH
      if (err/=0) then
        WRITE(u6,"('cannot allocate memory in preclinsol')")
        PRINT *,'cannot allocate memory in preclinsol. domain ',
     >                       name(1:idx(proc))
      endif
      
!$OMP BARRIER
!$OMP FLUSH

      CALL collecterr
      if (err/=0) return
      
#if 0
      IF (proc==1.and.ISWEEPa.EQ.1.and.ISTAT1.LE.251.and.
     >    251.LT.ISTAT1+l1b*ISTEP) THEN
            CALL makecrs(PP, Q, R, S, T,
     >             UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
      END IF  
#endif      
C
C  restrtgmr:  restart parameter for GMRES.
C  lmax:  restart for BICGSTABl.
C  maxliniter: the maximum number of GMRES iterations to be performed.
C  iabs:       determines stopping criterion for GMRES.
C  ioutgmr:    determines how much output is printed, used in GMRES.
C  nonzerox:   determines whether or not initial quess is zero. If
C              this is the case, the number of matrix-vector multiplications
C              reduces by one.
C 
      IF (abs(relaxall).GT.1.D-10) 
     >    CALL underrel(b,x,PP, Q, R, S, T,
     >                  UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
C  ------------------------
C  COMPUTE INITIAL RESIDUAL
C  ------------------------
C
      ioutgmr=0
      nummatvec=nummatvec+1
      CALL matvec(work,x,va,Ycomm, PP, Q, R, S, T,
     >            UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
      if (err/=0) return

      b(1:Nt)=b(1:Nt)-work(1:Nt)

#if 0
      IF (proc==1.and.ISWEEPa.EQ.1.and.ISTAT1.LE.71.and.
     >    71.LT.ISTAT1+l1b*ISTEP) THEN
        OPEN (47, FILE = 'righthandside.dat')
        REWIND (47)
        write(47,*) Nt
        DO i=1,Nt
          WRITE (47,*)  b(i)
        ENDDO
        CLOSE (47)
      END IF  
#endif

      CALL detmax(b,x,PP, Q, R, S, T,
     >            UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
      CALL collecterr
      if (err/=0) return
      
      IF (ISWEEPa==1.and.skipstation.lt.1.d-8) THEN
         CALL res2tec(b,PP, Q, R, S, T,
     >             UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
      END IF
      R2normloc=dnrm2own(Nt,b,1)
!      write(u6,*) 'Before R2normloc b',R2normloc
      CALL precon5p(b, PP, Q, R, S, T,
     >              e1,e5,e2,e4,
     >              UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
      CALL collecterr
      if (err/=0) return
      R2normloc=dnrm2own(Nt,b,1)
!      write(u6,*) 'After R2normloc b',R2normloc

      IF (rightprecon) THEN
      ELSE
C
C     Determine preconditioned solution.
C
        CALL solve5p(b,va,Ycomm, PP, Q, R, S, T,
     >               e1,e5,e2,e4,
     >               UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
      END IF
      if (knoll) then
        x(1:Nt)=b(1:Nt)
      else
        x(1:Nt)=zero
      endif
      nonzerox=knoll
C     R2normprecloc=dnrm2own(Nt,b,1)
      R2normprecloc=ZERO
      DO i=1,Nt
        R2normprecloc=R2normprecloc+b(i)*b(i)
      ENDDO
      R2normprecloc=DSQRT(R2normprecloc)
      Fraction=R2normprecloc/R2normloc

      scalararray(proc)=Fraction
!$OMP BARRIER
!$OMP FLUSH
      DO nb1=minproc,maxproc
        Fraction=MAX(Fraction,scalararray(nb1))
      ENDDO
!$OMP BARRIER
!$OMP FLUSH

      IF (rightprecon) THEN
      ELSE
        IF (Fraction.GT.1.D2) THEN
          WRITE(u6,"('!!! Fraction norm preconditioned residual and unpr
     &econditioned residual',1PE24.10)") Fraction
          WRITE(u6,*) 'ISTAT1 ',ISTAT1
        END IF
      END IF

      max2norm=MAX(max2norm,R2normprecloc)
      scalararray(proc)=max2norm
!$OMP BARRIER
!$OMP FLUSH
      DO nb1=minproc,maxproc
        max2norm=MAX(max2norm,scalararray(nb1))
      ENDDO
!$OMP BARRIER
!$OMP FLUSH

      R2normprec=R2normprec+R2normprecloc**2
C 
C  -------------------------
C  COMPUTE CORRECTION VECTOR
C  -------------------------
C
      liniter = maxliniter
      restrt = restrtgmr 
      IF (stopcritgmr.EQ.'abs') THEN
        tol1= tolgmr*max(max2normold,max2norm)
        iabs=1
      ELSE
C
C Acceptance criterion for GMRES includes estimate of the preconditioner.
C The constants below have been chosen in such a way the same input-files as in
C version_2.d00 can be used.
C
       tol1= 1.5D1*tolgmr/MAX(ONE,Fraction)
!june2012        tol1=tolgmr
        iabs=0
      END IF
      IF (MOD(ISWEEP,NUMOUTP).EQ.0) ioutgmr=0 
      
C
C       APPLY GMRES.
C
      CALL gmres(Nt,restrt,b,x,work,tol1,iabs,liniter,
     &           minlingmr,matvec,ope,ioutgmr,nonzerox,va,Ycomm, 
     >           PP, Q, R, S, T,e1,e5,e2,e4,
     >           UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb,NY,NZ,L1b,
     >           rightprecon,u6)
      tol2=tol1
      IF (stopcritgmr.EQ.'abs') tol1=tol1/R2normprecloc
      IF (MOD(ISWEEP,NUMOUTP).EQ.0) 
     &    WRITE(u6,"('    norm reduction in ',I4,' steps of gmres ',
     &                   1PE24.10)") liniter,tol1 
      nummatvec=nummatvec+liniter
        
      IF (rightprecon) THEN
C
C     Determine preconditioned solution.
C
        R2normloc=dnrm2own(Nt,x,1)
        CALL solve5p(x,va,Ycomm, PP, Q, R, S, T,
     >               e1,e5,e2,e4,
     >               UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
        R2normprecloc=dnrm2own(Nt,x,1)
        IF (MOD(ISWEEP,NUMOUTP).EQ.0) THEN
          IF (R2normprecloc.GT.1.D2*R2normloc) THEN
            WRITE(u6,"('!!! Fraction norm preconditioned solution and un
     &preconditioned solution',1PE24.10)") R2normprecloc/R2normloc 
          END IF
        END IF
      END IF

      if (nonzerox) nummatvec=nummatvec+1
      
      deallocate(work,e1,e5,e2,e4)

!      write(u6,*) 'tol2 ',tol2

!      IF (tol2.GT.1.d0) THEN
!        x(1:Nt)=ZERO
!        WRITE(u6,*) 'ISTAT1 ',ISTAT1,
!     >       ' Correction vector cannot be computed.'
!      END IF

      RETURN
      END

      REAL*8 FUNCTION ddotp(nn,xx,yy)
C
C Computes innerproduct of two vectors that are distributed
C over several processors.
C
      IMPLICIT NONE
      include "dataparal.incl"
      INTEGER  nn,nbt,ii
      REAL*8 xx(nn),yy(nn),dotpart,ddot,sumdot
C
      dotpart=0.D0
      DO ii=1,nn
        dotpart=dotpart+xx(ii)*yy(ii)
      ENDDO
C
C    Send innerprocduct of local vector to the other process.
C    Fetch innerprocduct of remote vector from the other process.
C
      scalararray(proc)=dotpart
!$OMP BARRIER
!$OMP FLUSH
      sumdot=0.D0
      DO nbt=minproc,maxproc
        sumdot=sumdot+scalararray(nbt)
      ENDDO
!$OMP BARRIER
!$OMP FLUSH
C
C     Compute innerproduct of the whole vector.
C
      ddotp=sumdot
      return
      end

      subroutine ope(nn,xx,yy,matvec,va,Ycomm, 
     >               PP, Q, R, S, T,
     >               e1,e5,e2,e4,
     >               UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb,
     >               rightprecon)
C
C     Auke van der Ploeg, CWI, 28-1-98
C
C Delivers matrix-vector multiplication yy := L*xx 
C by means of `matvec`
C
      include 'dataparnas.incl'
      include 'datamat.incl'
      include 'datava.incl'
      include 'dataprec.incl'
C
      LOGICAL rightprecon
      INTEGER    nn
      REAL*8 xx(nn),yy(nn)
      REAL*8,allocatable,dimension(:):: yy2

      external matvec

      if (rightprecon) then   
        allocate(yy2(nn))
        yy2(1:nn)=xx(1:nn) 
        CALL solve5p(yy2,va,Ycomm, PP, Q, R, S, T,
     >               e1,e5,e2,e4,
     >               UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
        CALL matvec(yy,yy2,va,Ycomm, PP, Q, R, S, T,
     >              UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
        deallocate(yy2)
      else
        CALL matvec(yy,xx,va,Ycomm, PP, Q, R, S, T,
     >              UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
        CALL solve5p(yy,va,Ycomm, PP, Q, R, S, T,
     >               e1,e5,e2,e4,
     >               UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
      endif
      RETURN
      END
C
      subroutine matvec(y,x,va,Ycomm, PP, Q, R, S, T,
     >                  UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
C
C     Auke van der Ploeg, April 9 1998.
C
*          Contains the matrix-vector multiplication
*
*               y := Ax.
C
      include 'dataparnas.incl'
      include "dataparal.incl"
      INCLUDE 'datamat.incl'
      include 'datava.incl'
      INTEGER uu,vv,kop
      integer nbJ,jb,jbA,nbnd,nbK,kb,kbA,step1(6)
      REAL*8 x(N,1:NY,1:ubk,L1b), y(N,1:NY,1:ubk,L1b), t1 
C
      real*8 ddot,ff
      integer cnt,cntbew,KR
      LOGICAL bool

      t1=dsecnd()
      Ycomm(1:N,1:max(NY,ubk),1:L1b,3:6)=ZERO
      IF (maxproc.EQ.minproc) GO TO 100
C----------------------------------------------------
C    Communication across a j-is-constant line.
C----------------------------------------------------

      step1(1:N)=1
      step1(5:6)=KSTEP
      DO L1=1,L1b
      ISTAT=ISTAT1+ISTEP*(L1-1)
      va(1:N,0:NY+1,0)=ZERO
      va(1:N,0:NY+1,ubk+1)=ZERO
      va(1:N,0,0:ubk+1)=ZERO
      va(1:N,NY+1,0:ubk+1)=ZERO
      va(1:N,1:NY,1:ubk)=x(1:N,1:NY,1:ubk,L1)
!$OMP BARRIER
!$OMP FLUSH
      DO nbnd=3,6
        nbJ=commtableL1(L1,nbnd,1)
        IF (nbJ.GT.0) THEN
          trx(1:N*ubAop(nbnd,L1))=0.d0
          DO uu=0,1
            DO vv=0,1
              cnt=0
                DO KW=1,ubAop(nbnd,L1)
                  kop=1+(KW-1)*step1(nbnd)
                  jm=jindexnc(kop,L1,nbnd)-1
                  km=kindexnc(kop,L1,nbnd)-1
                  jj=jm+uu
                  kk=km+vv
                  trx(cnt+1:cnt+N)=trx(cnt+1:cnt+N)+
     >              weights(uu,vv,kop,L1,nbnd)*va(1:N,jj,kk)
                  cnt=cnt+N
                ENDDO
            ENDDO
          ENDDO
          trxshared(1:cnt,commtableL1(L1,nbnd,2),nbJ)=trx(1:cnt)
        END IF
      ENDDO
!$OMP BARRIER
!$OMP FLUSH
      DO nbnd=5,6
        nbJ=commtableL1(L1,nbnd,1)
        IF (nbJ.GT.0) THEN
          cnt=N*ubk
          trx(1:cnt)=trxshared(1:cnt,nbnd,proc)
          cnt=0
            DO k=1,ubk
              KR=1+(k-1)*KSTEP
              IF (nbnd.NE.5.or.IBCJ1(KR,ISTAT).EQ.10.
     >                      or.IBCJ1(KR,ISTAT).EQ.5) THEN
                Ycomm(1:N,k,L1,nbnd)=trx(cnt+1:cnt+N)
              END IF
              cnt=cnt+N
            ENDDO
        END IF
      ENDDO
!$OMP BARRIER
!$OMP FLUSH
      DO nbnd=3,4
        nbk=commtableL1(L1,nbnd,1)
        IF (nbk.GT.0) THEN
           DO j=1,NY
             Ycomm(1:N,j,L1,nbnd)=trxshared(N*j-N+1:N*j,nbnd,proc)
           ENDDO
        END IF
      ENDDO
!$OMP BARRIER
!$OMP FLUSH
C
C  Loop over L1:
C
      ENDDO
C
C     Matrix-vector multiplication y=Ax is performed:
C
 100  DO L1=1,L1b
        ISTAT=ISTAT1+ISTEP*(L1-1)
        va(1:N,1:NY,1:ubk)=x(1:N,1:NY,1:ubk,L1)
        va(1:N,objAL1(L1)-1,1:ubk)=Ycomm(1:N,1:ubk,L1,5)
        va(1:N,NY+1,1:ubk)=Ycomm(1:N,1:ubk,L1,6)
        va(1:N,1:NY,obkAL1(L1)-1)=Ycomm(1:N,1:NY,L1,3)
        va(1:N,1:NY,ubk+1)=Ycomm(1:N,1:NY,L1,4)
        DO k=1,ubk
          DO j=1,NY
            y(1:N,j,k,L1)=
     &                Q (1:N,1,j,k,L1)*va(1,j,k)
     &               +Q (1:N,2,j,k,L1)*va(2,j,k)
     &               +Q (1:N,3,j,k,L1)*va(3,j,k)
     &               +Q (1:N,4,j,k,L1)*va(4,j,k)
     &               +S (1:N,1,j,k,L1)*va(1,j,k-1)
     &               +S (1:N,2,j,k,L1)*va(2,j,k-1)
     &               +S (1:N,3,j,k,L1)*va(3,j,k-1)
     &               +S (1:N,4,j,k,L1)*va(4,j,k-1)
     &               +T (1:N,1,j,k,L1)*va(1,j,k+1)
     &               +T (1:N,2,j,k,L1)*va(2,j,k+1)
     &               +T (1:N,3,j,k,L1)*va(3,j,k+1)
     &               +T (1:N,4,j,k,L1)*va(4,j,k+1)
     &               +PP(1:N,1,j,k,L1)*va(1,j-1,k)
     &               +PP(1:N,2,j,k,L1)*va(2,j-1,k)
     &               +PP(1:N,3,j,k,L1)*va(3,j-1,k)
     &               +PP(1:N,4,j,k,L1)*va(4,j-1,k)
     &               +R (1:N,1,j,k,L1)*va(1,j+1,k)
     &               +R (1:N,2,j,k,L1)*va(2,j+1,k)
     &               +R (1:N,3,j,k,L1)*va(3,j+1,k)
     &               +R (1:N,4,j,k,L1)*va(4,j+1,k)     
C    &          MATMUL(Q (1:N,1:N,j,k,L1),va(1:N,j,k))
C    &         +MATMUL(S (1:N,1:N,j,k,L1),va(1:N,j,k-1))
C    &         +MATMUL(T (1:N,1:N,j,k,L1),va(1:N,j,k+1))
C    &         +MATMUL(PP(1:N,1:N,j,k,L1),va(1:N,j-1,k))
C    &         +MATMUL(R (1:N,1:N,j,k,L1),va(1:N,j+1,k))
          ENDDO
        ENDDO
        IF (L1.EQ.2) THEN
          DO k=1,ubk
            DO j=1,NY
              y(1:N,j,k,2)=y(1:N,j,k,2)
     &                     +UPSTRa(1:N,1,j,k,2)*x(1,j,k,1)
     &                     +UPSTRa(1:N,2,j,k,2)*x(2,j,k,1)
     &                     +UPSTRa(1:N,3,j,k,2)*x(3,j,k,1)
     &                     +UPSTRa(1:N,4,j,k,2)*x(4,j,k,1)
C    &         +MATMUL(UPSTRa(1:N,1:N,j,k,2),x(1:N,j,k,1))
            ENDDO
          ENDDO
        END IF
        IF (L1.GT.2) THEN
          DO k=1,ubk
            DO j=1,NY
              y(1:N,j,k,L1)=y(1:N,j,k,L1)
     &                      +UPSTRa(1:N,1,j,k,L1)*x(1,j,k,L1-1)
     &                      +UPSTRa(1:N,2,j,k,L1)*x(2,j,k,L1-1)
     &                      +UPSTRa(1:N,3,j,k,L1)*x(3,j,k,L1-1)
     &                      +UPSTRa(1:N,4,j,k,L1)*x(4,j,k,L1-1)
     &                      +UPSTRb(1:N,1,j,k,L1)*x(1,j,k,L1-2)
     &                      +UPSTRb(1:N,2,j,k,L1)*x(2,j,k,L1-2)
     &                      +UPSTRb(1:N,3,j,k,L1)*x(3,j,k,L1-2)
     &                      +UPSTRb(1:N,4,j,k,L1)*x(4,j,k,L1-2)
C    &         +MATMUL(UPSTRa(1:N,1:N,j,k,L1),x(1:N,j,k,L1-1))
C    &         +MATMUL(UPSTRb(1:N,1:N,j,k,L1),x(1:N,j,k,L1-2))
            ENDDO
          ENDDO
        ENDIF
        IF (L1+1.EQ.L1b) THEN
          DO k=1,ubk
            DO j=1,NY
              y(1:N,j,k,L1)=y(1:N,j,k,L1)
     &                      +DOWNSTRa(1:N,1,j,k,L1)*x(1,j,k,L1+1)
     &                      +DOWNSTRa(1:N,2,j,k,L1)*x(2,j,k,L1+1)
     &                      +DOWNSTRa(1:N,3,j,k,L1)*x(3,j,k,L1+1)
     &                      +DOWNSTRa(1:N,4,j,k,L1)*x(4,j,k,L1+1)
C    &         +MATMUL(DOWNSTRa(1:N,1:N,j,k,L1),x(1:N,j,k,L1+1))
            ENDDO
          ENDDO
        END IF
        IF (L1+2.LE.L1b) THEN
          DO k=1,ubk
            DO j=1,NY
              y(1:N,j,k,L1)=y(1:N,j,k,L1)
     &                      +DOWNSTRa(1:N,1,j,k,L1)*x(1,j,k,L1+1)
     &                      +DOWNSTRa(1:N,2,j,k,L1)*x(2,j,k,L1+1)
     &                      +DOWNSTRa(1:N,3,j,k,L1)*x(3,j,k,L1+1)
     &                      +DOWNSTRa(1:N,4,j,k,L1)*x(4,j,k,L1+1)
     &                      +DOWNSTRb(1:N,1,j,k,L1)*x(1,j,k,L1+2)
     &                      +DOWNSTRb(1:N,2,j,k,L1)*x(2,j,k,L1+2)
     &                      +DOWNSTRb(1:N,3,j,k,L1)*x(3,j,k,L1+2)
     &                      +DOWNSTRb(1:N,4,j,k,L1)*x(4,j,k,L1+2)
C    &         +MATMUL(DOWNSTRa(1:N,1:N,j,k,L1),x(1:N,j,k,L1+1))
C    &         +MATMUL(DOWNSTRb(1:N,1:N,j,k,L1),x(1:N,j,k,L1+2))
            ENDDO
          ENDDO
        END IF
        if (objAL1(L1).GT.1) y(1:N,1,1:ubk,L1)=x(1:N,1,1:ubk,L1)
        if (obkAL1(L1).GT.1) y(1:N,1:NY,1,L1)=x(1:N,1:NY,1,L1)
!        if (objAL1(L1).GT.1) y(1:N,1,1:ubk,L1)=va(1:N,1,1:ubk)
!        if (obkAL1(L1).GT.1) y(1:N,1:NY,1,L1)=va(1:N,1:NY,1)
      ENDDO
      
            
      ttotmv=ttotmv+dsecnd()-t1
      RETURN
      end
C
      subroutine underrel(yy,xx,PP, Q, R, S, T,
     >                    UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
      include 'dataparnas.incl'
      INCLUDE 'datamat.incl'
C
      REAL*8, DIMENSION(N,1:NY,1:ubk,L1b):: yy,xx

      IF (relaxall.GT.zero) THEN
        DO L1=1,L1b
          I=ISTAT1+(L1-1)*ISTEP
          DO k=obkAL1(L1)+1,ubk-1
            DO j=objAL1(L1)+1,NY-1
              yy(1:N,j,k,L1)=yy(1:N,j,k,L1)+relaxall*
     >         (Q(1:N,1,j,k,L1)*xx(1,j,k,L1)+
     >          Q(1:N,2,j,k,L1)*xx(2,j,k,L1)+
     >          Q(1:N,3,j,k,L1)*xx(3,j,k,L1)+
     >          Q(1:N,4,j,k,L1)*xx(4,j,k,L1))
              Q(1:N,1:N,j,k,L1)=(ONE+relaxall)*Q(1:N,1:N,j,k,L1)
            ENDDO
          ENDDO
        ENDDO
      ELSE
!
! No implicit underrelaxation for mass equation or pressure
!
        DO L1=1,L1b
          I=ISTAT1+(L1-1)*ISTEP
          DO k=obkAL1(L1)+1,ubk-1
            DO j=objAL1(L1)+1,NY-1
              DO ii=1,4
                IF (ii.ne.2) THEN
                  yy(ii,j,k,L1)=yy(ii,j,k,L1)-relaxall*
     >             (Q(ii,1,j,k,L1)*xx(1,j,k,L1)+
     >              Q(ii,2,j,k,L1)*xx(2,j,k,L1)+
     >              Q(ii,3,j,k,L1)*xx(3,j,k,L1))
                END IF
              ENDDO
              Q(1,1:3,j,k,L1)=(ONE-relaxall)*Q(1,1:3,j,k,L1)
              Q(3,1:3,j,k,L1)=(ONE-relaxall)*Q(3,1:3,j,k,L1)
              Q(4,1:3,j,k,L1)=(ONE-relaxall)*Q(4,1:3,j,k,L1)
           ENDDO
          ENDDO
        ENDDO
!        DO L1=1,L1b
!          I=ISTAT1+(L1-1)*ISTEP
!          DO k=obkAL1(L1)+1,ubk-1
!            DO j=objAL1(L1)+1,NY-1
!              DO ii=2,2
!                  yy(ii,j,k,L1)=yy(ii,j,k,L1)-relaxall*
!     >             (Q(ii,1,j,k,L1)*xx(1,j,k,L1)+
!     >              Q(ii,2,j,k,L1)*xx(2,j,k,L1)+
!     >              Q(ii,3,j,k,L1)*xx(3,j,k,L1))
!              ENDDO
!              Q(2,1:3,j,k,L1)=(ONE-relaxall)*Q(2,1:3,j,k,L1)
!           ENDDO
!          ENDDO
!        ENDDO
      END IF
      RETURN
      END

      subroutine res2tec(yy,PP, Q, R, S, T,
     >                  UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
C
C     determine maximum residual and location.
C     Residual obtained after scaling
C     sum(abs(Q(ii,1..4))) to ONE.
C
      include 'dataparnas.incl'
      include 'dataparal.incl'
      INCLUDE 'datamat.incl'
C
      REAL*8, DIMENSION(N,1:NY,1:ubk,L1b):: yy
      REAL*8 ff(4),ff1
      INTEGER i1,tecu

      tecu=u21      
      DO L1=1,L1b
        I=ISTAT1+(L1-1)*ISTEP
        DO k=1,ubk
          DO j=1,NY
            DO i1=1,N
              ff1=abs(Q(i1,1,j,k,L1))+
     >            abs(Q(i1,2,j,k,L1))+
     >            abs(Q(i1,3,j,k,L1))+
     >            abs(Q(i1,4,j,k,L1))
              ff(i1)=4.D0*yy(i1,j,k,L1)/ff1
            ENDDO
            write(tecu,'(1P7E20.12)') GP(1:3,j,K,I),ff(1:4)
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END

      subroutine detmax(yy,xx,PP, Q, R, S, T,
     >                  UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
C
C     determine maximum residual and location.
C     Residual obtained after scaling
C     Q(ii,1..4) to ONE.
C
      include 'dataparnas.incl'
      include 'dataparal.incl'
      INCLUDE 'datamat.incl'
C
      REAL*8, DIMENSION(N,1:NY,1:ubk,L1b):: yy,xx
      REAL*8 maxF(1:N),ff,ff1
      INTEGER jloc(1:N),kloc(1:N),i1,ubktemp

      DO L1=1,L1b
        I=ISTAT1+(L1-1)*ISTEP
        maxF(1:N)=ZERO
        jloc(1:N)=0
        kloc(1:N)=0
        ubktemp=ubk
        IF (IBCKNZ(1,I)==10.and.
     >       (IBCK1(1,I)==6.or.IBCK1(1,I)==11)) ubktemp=ubk-1
        DO k=obkAL1(L1),ubk
          DO j=objAL1(L1),NY
            DO i1=1,N
              ff1=abs(Q(i1,1,j,k,L1))+
     >            abs(Q(i1,2,j,k,L1))+
     >            abs(Q(i1,3,j,k,L1))+
     >            abs(Q(i1,4,j,k,L1))
              IF (ff1.LT.1.D-20) THEN
                write(u6,*) 'Equation ',i1,' is zero.'
                write(u6,*) 'Point ijk ',I,j,k
                PRINT *,'Equation ',i1,' is zero. domain ',
     >                       name(1:idx(proc))
                PRINT *,'Point ijk ',I,j,k
                err=-1
                return
              END IF
              ff=4.D0*abs(yy(i1,j,k,L1)/ff1)
              R2norm=R2norm+ff**2
              IF (ff.gt.maxF(i1).and.k.LE.ubktemp) then
                maxF(i1)=ff
                jloc(i1)=j
                kloc(i1)=1+(k-1)*kstep
              END IF
              
!
! Extra implicit underrelaxation in case of divergence:
!              
!              if (ff.gt.2.d0) then
!                yy(i1,j,k,L1)=yy(i1,j,k,L1)+1.d0*
!     >           (Q(i1,1,j,k,L1)*xx(1,j,k,L1)+
!     >            Q(i1,2,j,k,L1)*xx(2,j,k,L1)+
!     >            Q(i1,3,j,k,L1)*xx(3,j,k,L1)+
!     >            Q(i1,4,j,k,L1)*xx(4,j,k,L1))
!                Q(i1,1:N,j,k,L1)=(ONE+ONE)*Q(i1,1:N,j,k,L1)
!              endif

              
              
!              IF (ff.gt.1.d20) THEN
!                write(u6,*),'i1 j k I Q yy ',i1,j,k,I,Q(i1,1:4,j,k,L1),
!     >                            yy(i1,j,k,L1)   
!              END IF
            ENDDO
          ENDDO
        ENDDO
        IF (MOD(ISWEEP,NUMOUTP).EQ.0) THEN
          WRITE(u6,8350) I,maxF(1),jloc(1),kloc(1),
     >                     maxF(2),jloc(2),kloc(2),
     >                     maxF(3),jloc(3),kloc(3),
     >                     maxF(4),jloc(4),kloc(4)
        END IF
        DO i1=1,N
          IF(MAXF(i1).GT.MAXR(i1)) THEN
            MAXR(i1)=maxF(i1)
            IPMAXR(i1)=I
            JPMAXR(i1)=jloc(i1)
            KPMAXR(i1)=kloc(i1)
          ENDIF
        ENDDO
      ENDDO

 8350 FORMAT(7X,'Max scaled res. I=',I3,': ksi-m. ',1PE10.3,
     >' at (J,K) = ('
     >,I3,',',I3,');  conti. ',1PE10.3,' at (J,K) = (',I3,',',I3,');'/,
     >7X,'                       zeta-m.',1PE10.3,' at (J,K) = (',I3,','
     >,I3,');  eta-m. ',1PE10.3,' at (J,K) = (',I3,',',I3,')')
      RETURN
      END

      subroutine solve5p(yy,va,Ycomm, PP, Q, R, S, T,
     >                   e1,e5,e2,e4,
     >                   UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
C
C     Auke van der Ploeg, August 7 1998.
C
C     This routine applies the preconditoner which  
C     must have been constructed in subroutine 'precon5p.' 
C     It solves the linear system
C         L U x = yy
C     in which  L  is a sparse block-lower-triangular matrix
C     and U  is a sparse block-upper-triangular matrix with identity
C     blocks on the diagonal. Put the solution back in yy.
C
C     Modified version, October 14 1999.
C
C     Performs block-triangular solves.
C     The upper-triangular blocks coming from the discretisation of
C     the pressure in xi-direction are also taken care of:
C
C
C     The block-pentadigonal matrix
C
C        d1 u1 v1
C        s2 d2 u2 v2
C        t3 s3 d3 u3 v3
C           t4 s4 d4 u4 v4
C             \  \  \  \  \
C    
C     is approximated by the block-lower-triangular matrix
C
C        d1 0
C        s2 d2 0 
C        t3 s3 d3 0 
C           t4 s4 d4 0  
C             \  \  \  \
C
C     multiplied by the block-upper triangular matrix
C
C        I inv(d1)u1 inv(d1)v1
C          I         inv(d2)u2 inv(d2)v2
C                    I         inv(d3)u3  inv(d3)v3
C                              I          inv(d4)u4  inv(d4)v4
C                               \          \         \          \
C
C     in which d1, d2, etc. are incomplete decompositions of the
C     blocks on the main diagonal.
C
      include 'dataparnas.incl'
      INCLUDE 'datamat.incl'
      INCLUDE 'datava.incl'
      INCLUDE 'dataprec.incl'
      REAL*8 yy(N,1:NY,1:ubk,L1b) 
      REAL*8 t1
C=======================================================================
C         Parameters:
C=======================================================================
C
C     yy      i/o On input: right hand side of the linear system.
C                 On output: solution vector
C
      va(1:N,0,0:ubk+1)=ZERO
      va(1:N,NY+1,0:ubk+1)=ZERO
      va(1:N,0:NY+1,0)=ZERO
      va(1:N,0:NY+1,ubk+1)=ZERO
      t1=dsecnd()
      DO L1=1,L1b
        IF (L1.EQ.2) THEN
          DO k=1,ubk
            DO j=1,NY
              yy(1:N,j,k,2)=yy(1:N,j,k,2)
     &                    -UPSTRa(1:4,1,j,k,2)*yy(1,j,k,1)
     &                    -UPSTRa(1:4,2,j,k,2)*yy(2,j,k,1)
     &                    -UPSTRa(1:4,3,j,k,2)*yy(3,j,k,1)
     &                    -UPSTRa(1:4,4,j,k,2)*yy(4,j,k,1)
C    &           -MATMUL(UPSTRa(1:N,1:N,j,k,2),yy(1:N,j,k,1))
            ENDDO
          ENDDO
        END IF
        IF (L1.GT.2) THEN
          DO k=1,ubk
            DO j=1,NY
              yy(1:N,j,k,L1)=yy(1:N,j,k,L1)
     &                     -UPSTRa(1:4,1,j,k,L1)*yy(1,j,k,L1-1)
     &                     -UPSTRa(1:4,2,j,k,L1)*yy(2,j,k,L1-1)
     &                     -UPSTRa(1:4,3,j,k,L1)*yy(3,j,k,L1-1)
     &                     -UPSTRa(1:4,4,j,k,L1)*yy(4,j,k,L1-1)
     &                     -UPSTRb(1:4,1,j,k,L1)*yy(1,j,k,L1-2)
     &                     -UPSTRb(1:4,2,j,k,L1)*yy(2,j,k,L1-2)
     &                     -UPSTRb(1:4,3,j,k,L1)*yy(3,j,k,L1-2)
     &                     -UPSTRb(1:4,4,j,k,L1)*yy(4,j,k,L1-2)
C    &           -MATMUL(UPSTRa(1:N,1:N,j,k,L1),yy(1:N,j,k,L1-1))
C    &           -MATMUL(UPSTRb(1:N,1:N,j,k,L1),yy(1:N,j,k,L1-2))
            ENDDO
          ENDDO
        END IF
C
C      va:=(I-L)^{-1}yy
C
      DO k=1,ubk
        DO j=1,NY
           va(1:N,j,k)=yy(1:N,j,k,L1) 
     &              - e2(1:4,1,j,k,0,L1)*va(1,j-1,k)
     &              - e2(1:4,2,j,k,0,L1)*va(2,j-1,k)
     &              - e2(1:4,3,j,k,0,L1)*va(3,j-1,k)
     &              - e2(1:4,4,j,k,0,L1)*va(4,j-1,k)
     &              - e1(1:4,1,j,k,1,L1)*va(1,j+1,k-1)
     &              - e1(1:4,2,j,k,1,L1)*va(2,j+1,k-1)
     &              - e1(1:4,3,j,k,1,L1)*va(3,j+1,k-1)
     &              - e1(1:4,4,j,k,1,L1)*va(4,j+1,k-1)
     &              - S (1:4,1,j,k,L1)*va(1,j  ,k-1)
     &              - S (1:4,2,j,k,L1)*va(2,j  ,k-1)
     &              - S (1:4,3,j,k,L1)*va(3,j  ,k-1)
     &              - S (1:4,4,j,k,L1)*va(4,j  ,k-1)
C    &           -MATMUL(e2(1:N,1:N,j,k,0,L1),va(1:N,j-1,k))
C    &           -MATMUL(e1(1:N,1:N,j,k,1,L1),va(1:N,j+1,k-1))
C    &           -MATMUL(S (1:N,1:N,j,k,L1),  va(1:N,j  ,k-1))
        ENDDO
      ENDDO
      DO k=ubk,1,-1
        DO j=NY,1,-1
          va(1:N,j,k) = va(1:N,j,k)
     &              - e4(1:4,1,j,k,0,L1)*va(1,j+1,k)
     &              - e4(1:4,2,j,k,0,L1)*va(2,j+1,k)
     &              - e4(1:4,3,j,k,0,L1)*va(3,j+1,k)
     &              - e4(1:4,4,j,k,0,L1)*va(4,j+1,k)
     &              - e5(1:4,1,j,k,1,L1)*va(1,j-1,k+1)
     &              - e5(1:4,2,j,k,1,L1)*va(2,j-1,k+1)
     &              - e5(1:4,3,j,k,1,L1)*va(3,j-1,k+1)
     &              - e5(1:4,4,j,k,1,L1)*va(4,j-1,k+1)
     &              - T (1:4,1,j,k,L1)*va(1,j  ,k+1)
     &              - T (1:4,2,j,k,L1)*va(2,j  ,k+1)
     &              - T (1:4,3,j,k,L1)*va(3,j  ,k+1)
     &              - T (1:4,4,j,k,L1)*va(4,j  ,k+1)
C    &           -MATMUL(e4(1:N,1:N,j,k,0,L1),va(1:N,j+1,k))
C    &           -MATMUL(e5(1:N,1:N,j,k,1,L1),va(1:N,j-1,k+1))
C    &           -MATMUL(T (1:N,1:N,j,k,L1),  va(1:N,j  ,k+1))
        ENDDO
      ENDDO
      yy(1:N,1:NY,1:ubk,L1)=va(1:N,1:NY,1:ubk)
      ENDDO
      ttotsol=ttotsol+dsecnd()-t1
      RETURN
C
      END
C
      subroutine precon5p(ZZ1, PP, Q, R, S, T,
     >                    e1,e5,e2,e4,
     >                    UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb)
C
C     Auke van der Ploeg, August 7 1998.
C
C     This routine performs the preprocessing necessary for a 
C     preconditioner that can be used in PARNASSOS.
C     The matrices in PARNASSOS have a block-pentadigonal structure
C
C        d1 u1 v1
C        s2 d2 u2 v2
C        t3 s3 d3 u3 v3
C           t4 s4 d4 u4 v4
C             \  \  \  \  \
C
C     As preconditioner for GMRES we use the block-lower triangular matrix
C
C        d1 0
C        s2 d2 0
C        t3 s3 d3 0
C           t4 s4 d4 0
C             \  \  \  \
C
C     multiplied by the block-upper triangular matrix
C
C        I inv(d1)u1 inv(d1)v1
C          I         inv(d2)u2 inv(d2)v2
C                    I         inv(d3)u3  inv(d3)v3
C                              I          inv(d4)u4  inv(d4)v4
C                               \          \         \          \
C
C     in which d1, d2, etc. are incomplete decompositions of the
C     blocks on the main diagonal. These incomplete decompositions are
C     constructed in this routine.
C
C     Note that this routine can be parallellized, since all
C     incomplete decompositions of di can be computed simultaneously. 
C
C     The incomplete LU-decomposition of d1, d2 etc. are performed on
C     a 4*4-block level. 
*
*          See, for example, page 19 of the Phd thesis:
*          'Preconditioning for sparse matrices with applications.'
*          1994, University of Groningen.
*
      include 'dataparnas.incl'
      include 'dataparal.incl'
      INCLUDE 'datamat.incl'
      INCLUDE 'dataprec.incl'
      REAL*8 ZZ1(N,1:NY,1:ubk,L1b)
      integer cnt
C
      REAL*8 t1
C=======================================================================
C         Local variables:
C=======================================================================
      REAL*8 ff,b(N,N),inv(N,N),e11(N,N),e12(N,N),e20(N,N),relaxGus
      REAL*8 bupa(N,N),bupb(N,N),bdowna(N,N),bdownb(N,N),
     &       bQ(N,N),bP(N,N),bR(N,N),bS(N,N),bT(N,N),be4(N,N),be5(N,N)
      INTEGER num,jp2
C
      relaxGus=one/(one+relaxprecon)
      t1=dsecnd()
      DO L1=1,L1b
C
C  Initialization of e4 and e5 at k=1:
C
        e4(1:N,1:N,1:NY,1,0,L1)=ZERO
        e5(1:N,1:N,1:NY,1,1:1,L1)=ZERO
C
C  Initialization of e4 and e5 at j=1:
C
        e4(1:N,1:N,1,1:ubk,0,L1)=ZERO
        e5(1:N,1:N,1,1:ubk,1:1,L1)=ZERO
        DO k=1,ubk
          km=max(k-1,1)
          DO j=1,NY
            jm=max(j-1,1)
            jp=min(j+1,NY)
            jp2=min(j+2,NY)
            bupa(1:N,1:N)=UPSTRa(1:N,1:N,j,k,L1)
            bdowna(1:N,1:N)=DOWNSTRa(1:N,1:N,j,k,L1)
            bupb(1:N,1:N)=UPSTRb(1:N,1:N,j,k,L1)
            bdownb(1:N,1:N)=DOWNSTRb(1:N,1:N,j,k,L1)
            bQ(1:N,1:N)=Q(1:N,1:N,j,k,L1)
            bP(1:N,1:N)=PP(1:N,1:N,j,k,L1)
            bR(1:N,1:N)=R(1:N,1:N,j,k,L1)
            bS(1:N,1:N)=S(1:N,1:N,j,k,L1)
            bT(1:N,1:N)=T(1:N,1:N,j,k,L1)
            DO kk=1,N
              e11(1:N,kk)= - bS(1:N,1)*e4(1,kk,j,km,0,L1)
     &                     - bS(1:N,2)*e4(2,kk,j,km,0,L1)
     &                     - bS(1:N,3)*e4(3,kk,j,km,0,L1)
     &                     - bS(1:N,4)*e4(4,kk,j,km,0,L1)
              e20(1:N,kk)=bP(1:N,kk)
     &                - bS(1:N,1)*e5(1,kk,j,km,1,L1)
     &                - bS(1:N,2)*e5(2,kk,j,km,1,L1)
     &                - bS(1:N,3)*e5(3,kk,j,km,1,L1)
     &                - bS(1:N,4)*e5(4,kk,j,km,1,L1)
              b(1:N,kk)=bQ(1:N,kk)
     &                - bS(1:N,1)*T (1,kk,j,km,L1)
     &                - bS(1:N,2)*T (2,kk,j,km,L1)
     &                - bS(1:N,3)*T (3,kk,j,km,L1)
     &                - bS(1:N,4)*T (4,kk,j,km,L1)
            ENDDO
            IF (k.EQ.1.or.j.EQ.NY) e11(1:N,1:N)=ZERO
            DO kk=1,N
              e12(1:N,kk)= - e11(1:N,1)*e4(1,kk,jp,km,0,L1)
     &                     - e11(1:N,2)*e4(2,kk,jp,km,0,L1)
     &                     - e11(1:N,3)*e4(3,kk,jp,km,0,L1)
     &                     - e11(1:N,4)*e4(4,kk,jp,km,0,L1)
            ENDDO

            IF (k.EQ.1.or.j.GE.NY-1) e12(1:N,1:N)=ZERO
            e11(1:N,1:N)= e11(1:N,1:N)+e12(1:N,1:N)
            IF (j.EQ.1) e20(1:N,1:N)=ZERO

            DO kk=1,N
              b(1:N,kk)=b(1:N,kk)
     &                - e11(1:N,1)*e5(1,kk,jp,km,1,L1)
     &                - e11(1:N,2)*e5(2,kk,jp,km,1,L1)
     &                - e11(1:N,3)*e5(3,kk,jp,km,1,L1)
     &                - e11(1:N,4)*e5(4,kk,jp,km,1,L1)
     &                - e20(1:N,1)*e4(1,kk,jm,k,0,L1)
     &                - e20(1:N,2)*e4(2,kk,jm,k,0,L1)
     &                - e20(1:N,3)*e4(3,kk,jm,k,0,L1)
     &                - e20(1:N,4)*e4(4,kk,jm,k,0,L1)
              be4(1:N,kk)=bR(1:N,kk)
     &                - e11(1:N,1)*T (1,kk,jp,km,L1)
     &                - e11(1:N,2)*T (2,kk,jp,km,L1)
     &                - e11(1:N,3)*T (3,kk,jp,km,L1)
     &                - e11(1:N,4)*T (4,kk,jp,km,L1)
              be5(1:N,kk)=
     &         - e20(1:N,1)*(e5(1,kk,jm,k,1,L1)+T (1,kk,jm,k,L1))
     &         - e20(1:N,2)*(e5(2,kk,jm,k,1,L1)+T (2,kk,jm,k,L1))
     &         - e20(1:N,3)*(e5(3,kk,jm,k,1,L1)+T (3,kk,jm,k,L1))
     &         - e20(1:N,4)*(e5(4,kk,jm,k,1,L1)+T (4,kk,jm,k,L1))
            ENDDO
            IF (j.EQ.NY) be4=ZERO
            IF (k.EQ.ubk.or.j.EQ.1) be5=ZERO
C
C   Cramer's rule (not efficient for N much larger than 4). 
C
          inv(1,1)=b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))
     &            -b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))
     &            +b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))
          inv(2,1)=b(2,3)*(b(3,4)*b(4,1)-b(3,1)*b(4,4))
     &            -b(2,4)*(b(3,3)*b(4,1)-b(3,1)*b(4,3))
     &            +b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))
          inv(3,1)=b(2,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
     &            -b(2,1)*(b(3,4)*b(4,2)-b(3,2)*b(4,4))
     &            +b(2,2)*(b(3,4)*b(4,1)-b(3,1)*b(4,4))
          inv(4,1)=b(2,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))
     &            -b(2,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))
     &            +b(2,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
          inv(1,2)=b(3,2)*(b(4,3)*b(1,4)-b(4,4)*b(1,3))
     &            -b(3,3)*(b(4,2)*b(1,4)-b(4,4)*b(1,2))
     &            +b(3,4)*(b(4,2)*b(1,3)-b(4,3)*b(1,2))
          inv(2,2)=b(3,3)*(b(4,4)*b(1,1)-b(4,1)*b(1,4))
     &            -b(3,4)*(b(4,3)*b(1,1)-b(4,1)*b(1,3))
     &            +b(3,1)*(b(4,3)*b(1,4)-b(4,4)*b(1,3))
          inv(3,2)=b(3,4)*(b(4,1)*b(1,2)-b(4,2)*b(1,1))
     &            -b(3,1)*(b(4,4)*b(1,2)-b(4,2)*b(1,4))
     &            +b(3,2)*(b(4,4)*b(1,1)-b(4,1)*b(1,4))
          inv(4,2)=b(3,1)*(b(4,2)*b(1,3)-b(4,3)*b(1,2))
     &            -b(3,2)*(b(4,1)*b(1,3)-b(4,3)*b(1,1))
     &            +b(3,3)*(b(4,1)*b(1,2)-b(4,2)*b(1,1))
          inv(1,3)=b(4,2)*(b(1,3)*b(2,4)-b(1,4)*b(2,3))
     &            -b(4,3)*(b(1,2)*b(2,4)-b(1,4)*b(2,2))
     &            +b(4,4)*(b(1,2)*b(2,3)-b(1,3)*b(2,2))
          inv(2,3)=b(4,3)*(b(1,4)*b(2,1)-b(1,1)*b(2,4))
     &            -b(4,4)*(b(1,3)*b(2,1)-b(1,1)*b(2,3))
     &            +b(4,1)*(b(1,3)*b(2,4)-b(1,4)*b(2,3))
          inv(3,3)=b(4,4)*(b(1,1)*b(2,2)-b(1,2)*b(2,1))
     &            -b(4,1)*(b(1,4)*b(2,2)-b(1,2)*b(2,4))
     &            +b(4,2)*(b(1,4)*b(2,1)-b(1,1)*b(2,4))
          inv(4,3)=b(4,1)*(b(1,2)*b(2,3)-b(1,3)*b(2,2))
     &            -b(4,2)*(b(1,1)*b(2,3)-b(1,3)*b(2,1))
     &            +b(4,3)*(b(1,1)*b(2,2)-b(1,2)*b(2,1))
          inv(1,4)=b(1,2)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))
     &            -b(1,3)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))
     &            +b(1,4)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))
          inv(2,4)=b(1,3)*(b(2,4)*b(3,1)-b(2,1)*b(3,4))
     &            -b(1,4)*(b(2,3)*b(3,1)-b(2,1)*b(3,3))
     &            +b(1,1)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))
          inv(3,4)=b(1,4)*(b(2,1)*b(3,2)-b(2,2)*b(3,1))
     &            -b(1,1)*(b(2,4)*b(3,2)-b(2,2)*b(3,4))
     &            +b(1,2)*(b(2,4)*b(3,1)-b(2,1)*b(3,4))
          inv(4,4)=b(1,1)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))
     &            -b(1,2)*(b(2,1)*b(3,3)-b(2,3)*b(3,1))
     &            +b(1,3)*(b(2,1)*b(3,2)-b(2,2)*b(3,1))
          ff =  (b(1,1)*inv(1,1) - b(1,2)*inv(2,1)
     &          +b(1,3)*inv(3,1) - b(1,4)*inv(4,1))
           
          IF (ABS(ff).LT.1.D-80) THEN
            WRITE(u6,"('ff in precon5p j k L1',1PE24.10,3I4)") ff,j,k,L1
            WRITE(u6,"('b11 ',1PE24.10,' b21 ',1PE24.10)") b(1,1),b(2,1)
            WRITE(u6,"('b12 ',1PE24.10,' b22 ',1PE24.10)") b(1,2),b(2,2)
            WRITE(u6,"('b13 ',1PE24.10,' b23 ',1PE24.10)") b(1,3),b(2,3)
            WRITE(u6,"('b14 ',1PE24.10,' b24 ',1PE24.10)") b(1,4),b(2,4)
            WRITE(u6,"('b31 ',1PE24.10,' b41 ',1PE24.10)") b(3,1),b(4,1)
            WRITE(u6,"('b32 ',1PE24.10,' b42 ',1PE24.10)") b(3,2),b(4,2)
            WRITE(u6,"('b33 ',1PE24.10,' b43 ',1PE24.10)") b(3,3),b(4,3)
            WRITE(u6,"('b34 ',1PE24.10,' b44 ',1PE24.10)") b(3,4),b(4,4)
            WRITE(u6,"(' ')")
            WRITE(u6,"('Q11 ',1PE24.10,' Q21 ',1PE24.10)") 
     >                     Q(1,1,j,k,L1),Q(2,1,j,k,L1) 
            WRITE(u6,"('Q12 ',1PE24.10,' Q22 ',1PE24.10)") 
     >                     Q(1,2,j,k,L1),Q(2,2,j,k,L1) 
            WRITE(u6,"('Q13 ',1PE24.10,' Q23 ',1PE24.10)") 
     >                     Q(1,3,j,k,L1),Q(2,3,j,k,L1) 
            WRITE(u6,"('Q14 ',1PE24.10,' Q24 ',1PE24.10)") 
     >                     Q(1,4,j,k,L1),Q(2,4,j,k,L1) 
            WRITE(u6,"('Q31 ',1PE24.10,' Q41 ',1PE24.10)") 
     >                     Q(3,1,j,k,L1),Q(4,1,j,k,L1) 
            WRITE(u6,"('Q32 ',1PE24.10,' Q42 ',1PE24.10)") 
     >                     Q(3,2,j,k,L1),Q(4,2,j,k,L1) 
            WRITE(u6,"('Q33 ',1PE24.10,' Q43 ',1PE24.10)") 
     >                     Q(3,3,j,k,L1),Q(4,3,j,k,L1) 
            WRITE(u6,"('Q34 ',1PE24.10,' Q44 ',1PE24.10)") 
     >                     Q(3,4,j,k,L1),Q(4,4,j,k,L1) 
            PRINT *,'Zero pivot block. domain ',name(1:idx(proc))
            err=-1
            return
          ELSE
            ff=ONE/ff
          END IF
          inv(1,1)= ff*inv(1,1)
          inv(2,1)=-ff*inv(2,1)
          inv(3,1)= ff*inv(3,1)
          inv(4,1)=-ff*inv(4,1)
          inv(1,2)=-ff*inv(1,2)
          inv(2,2)= ff*inv(2,2)
          inv(3,2)=-ff*inv(3,2)
          inv(4,2)= ff*inv(4,2)
          inv(1,3)= ff*inv(1,3)
          inv(2,3)=-ff*inv(2,3)
          inv(3,3)= ff*inv(3,3)
          inv(4,3)=-ff*inv(4,3)
          inv(1,4)=-ff*inv(1,4)
          inv(2,4)= ff*inv(2,4)
          inv(3,4)=-ff*inv(3,4)
          inv(4,4)= ff*inv(4,4)
C
C  To have a 'relaxed' form of incomplete LU-decomposition,
C  increase main diagonal slightly (controlled by relaxprecon).
C  This increases robustness, but is in general less efficient.
C
          inv(1:N,1:N)=inv(1:N,1:N)*relaxGus
C
C Scale A and right-hand side in such a way that
C diagonal of preconditioner becomes I.
C
          ZZ1(1:N,j,k,L1)=inv(1:N,1)*ZZ1(1,j,k,L1)
     $                   +inv(1:N,2)*ZZ1(2,j,k,L1)
     $                   +inv(1:N,3)*ZZ1(3,j,k,L1)
     $                   +inv(1:N,4)*ZZ1(4,j,k,L1)
          DO kk=1,4
             UPSTRa(1:N,kk,j,k,L1)=inv(1:N,1)*bupa(1,kk)
     $                            +inv(1:N,2)*bupa(2,kk)
     $                            +inv(1:N,3)*bupa(3,kk)
     $                            +inv(1:N,4)*bupa(4,kk)
             DOWNSTRa(1:N,kk,j,k,L1)=inv(1:N,1)*bdowna(1,kk)
     $                              +inv(1:N,2)*bdowna(2,kk)
     $                              +inv(1:N,3)*bdowna(3,kk)
     $                              +inv(1:N,4)*bdowna(4,kk)
             UPSTRb(1:N,kk,j,k,L1)=inv(1:N,1)*bupb(1,kk)
     $                            +inv(1:N,2)*bupb(2,kk)
     $                            +inv(1:N,3)*bupb(3,kk)
     $                            +inv(1:N,4)*bupb(4,kk)
             DOWNSTRb(1:N,kk,j,k,L1)=inv(1:N,1)*bdownb(1,kk)
     $                              +inv(1:N,2)*bdownb(2,kk)
     $                              +inv(1:N,3)*bdownb(3,kk)
     $                              +inv(1:N,4)*bdownb(4,kk)
             Q(1:N,kk,j,k,L1)=inv(1:N,1)*bQ(1,kk)
     $                       +inv(1:N,2)*bQ(2,kk)
     $                       +inv(1:N,3)*bQ(3,kk)
     $                       +inv(1:N,4)*bQ(4,kk)
             PP(1:N,kk,j,k,L1)=inv(1:N,1)*bP(1,kk)
     $                        +inv(1:N,2)*bP(2,kk)
     $                        +inv(1:N,3)*bP(3,kk)
     $                        +inv(1:N,4)*bP(4,kk)
             R(1:N,kk,j,k,L1)=inv(1:N,1)*bR(1,kk)
     $                       +inv(1:N,2)*bR(2,kk)
     $                       +inv(1:N,3)*bR(3,kk)
     $                       +inv(1:N,4)*bR(4,kk)
             S(1:N,kk,j,k,L1)=inv(1:N,1)*bS(1,kk)
     $                       +inv(1:N,2)*bS(2,kk)
     $                       +inv(1:N,3)*bS(3,kk)
     $                       +inv(1:N,4)*bS(4,kk)
             T(1:N,kk,j,k,L1)=inv(1:N,1)*bT(1,kk)
     $                       +inv(1:N,2)*bT(2,kk)
     $                       +inv(1:N,3)*bT(3,kk)
     $                       +inv(1:N,4)*bT(4,kk)
             e1(1:N,kk,j,k,1,L1)=inv(1:N,1)*e11(1,kk)
     $                          +inv(1:N,2)*e11(2,kk)
     $                          +inv(1:N,3)*e11(3,kk)
     $                          +inv(1:N,4)*e11(4,kk)
             e2(1:N,kk,j,k,0,L1)=inv(1:N,1)*e20(1,kk)
     $                          +inv(1:N,2)*e20(2,kk)
     $                          +inv(1:N,3)*e20(3,kk)
     $                          +inv(1:N,4)*e20(4,kk)
             e4(1:N,kk,j,k,0,L1)=inv(1:N,1)*be4(1,kk)
     $                          +inv(1:N,2)*be4(2,kk)
     $                          +inv(1:N,3)*be4(3,kk)
     $                          +inv(1:N,4)*be4(4,kk)
             e5(1:N,kk,j,k,1,L1)=inv(1:N,1)*be5(1,kk)
     $                          +inv(1:N,2)*be5(2,kk)
     $                          +inv(1:N,3)*be5(3,kk)
     $                          +inv(1:N,4)*be5(4,kk)
          ENDDO
        ENDDO 
      ENDDO
      ENDDO
      timeprec=timeprec+dsecnd()-t1
      RETURN
      END

      subroutine gmres(n,im,rhs,sol,ss,eps,iabs,maxits,minits,
     $                 matvec, ope, ioutgmr, nonzerox,va,Ycomm, 
     &                 PP, Q, R, Sm, Tm, e1,e5,e2,e4,
     >                 UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb,NY,NZ,L1b,
     >                 rightprecon,u6)
c*************************************************************
C This code was initially written by Youcef Saad
C then revised by Henk A. van der Vorst  
C and Mike Botchev (oct. 1996)
C ************************************************************ 
c gmres algorithm . simple version .  (may 23, 1985)
c parameter list:
c n      == size of problem
c im     == size of krylov subspace:  should not exceed 50 in this
c          version (can be reset in code. looking at comment below)
c rhs    == right hand side
c sol    == initial guess on input, approximate solution on output
c ss     == work space of size n times (maxdim+2)
c eps    == tolerance for stopping criterion. process is stopped
c           as soon as ( ||.|| is the euclidean norm):
c           || current residual||/||initial residual|| <= eps
c           on OUTPUT: actual achieved norm residual (if iabs.ne.0)
c           or achieved relative residual norm reduction
c iabs   == use absolute stopping criterion
c           if (iabs.ne.0) then process is stopped as soon as
c           || current residual || <= eps  
c maxits == maximum number of iterations allowed
c minits == minimum number of iterations allowed
c ioutgmr   == output unit number number for printing intermediate results
c           if (ioutgmr .le. 0) no statistics are printed.
c nonzerox  input parameter whichs determines whether or not the
c           initial guess is zero
c ----------------------------------------------------------------
c subroutines used =
c ope(n,x,y,matvec)  == matrix vector multiplication y <- A*x
c                       envoking subroutine matvec
c BLAS:
c ddotp        == dot product function
c dnrm2own       == euclidean norm function
c daxpy       == y <-- y+ax routine
c dscal       == x <-- ax routine
c dtrsv       == to solve linear system with a triangular matrix
c*************************************************************
c-------------------------------------------------------------
c arnoldi size should not exceed 50 in this version..
c to reset modify maxdim             ----------------
c-------------------------------------------------------------
      IMPLICIT NONE

      LOGICAL rightprecon
      INTEGER NZ,NY,L1b,u6
      INCLUDE 'dataprec.incl'
C     INCLUDE 'datamat.incl'
      REAL*8 PP(4,4,NY,NZ,L1b), Q(4,4,NY,NZ,L1b),
     >        R(4,4,NY,NZ,L1b), Sm(4,4,NY,NZ,L1b), 
     >       Tm(4,4,NY,NZ,L1b),
     >       UPSTRa(4,4,NY,NZ,L1b),     UPSTRb(4,4,NY,NZ,L1b),
     >       DOWNSTRa(4,4,NY,NZ,L1b), DOWNSTRb(4,4,NY,NZ,L1b), 
     >       va(4,0:NY+1,0:NZ+1), Ycomm(4,MAX(NY,NZ),L1b,3:6)
      REAL*8 one
      LOGICAL nonzerox
      integer maxdim, maxd1,minits,ii
      parameter (maxdim=150, maxd1=maxdim+1, one=1.D0)
      external matvec,ope
      integer i,i1,iabs,im,ioutgmr,its,j,k,k1,maxits,n    !+mb
      REAL*8 coeff,dabs,ddotp,dnrm2own,dsqrt,eps,eps1,epsmac, roinit,
     &                 gam,rhs(n),ro,sol(n),ss(n,im+1),t     !+mb  
      REAL*8 hh(maxd1,maxdim),c(maxd1),s(maxd1),rs(maxd1) !+mb
 
      data epsmac/1.d-64/
 
      if (im .gt. maxdim) then
         im = maxdim
         write (*,'(A,i2)') 'GMRES: dimension has been reduced to ',im
         write (*,'(A)') ' => reset MAXDIM if you want it to be more'
      end if
      its = 0
c-------------------------------------------------------------
c **  outer loop starts here..
c-------------- compute initial residual vector --------------
 10   continue
C
C           ss(1):=A*sol
C
         IF (nonzerox.or.its.GT.0) THEN
           call ope(n,sol,ss,matvec,va,Ycomm, PP, Q, R,
     >              Sm, Tm, e1,e5,e2,e4,
     >              UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb,
     >              rightprecon)
c          ss(1:n,1)=rhs-ss(1:n,1)
           do ii=1,n
             ss(ii,1)=rhs(ii)-ss(ii,1)
           enddo
         ELSE
c          ss(1:n,1)=rhs
           do ii=1,n
             ss(ii,1)=rhs(ii)
           enddo
         END IF
c-------------------------------------------------------------
         ro = ddotp(n, ss, ss)
         ro = dsqrt(ro)
         roinit = ro
         if (ro .eq. 0.d0) then
            eps = ro
            maxits = its 
            return
         end if
 
         if (its .eq. 0) then      
C           set eps1 for stopping criterion
            if (iabs .ne. 0) then       !+mb use absolute criterion
               eps1=eps
               if (ro .le. eps1.and.minits.eq.0) then !+mb quit if accurate enough
                  eps = ro
                  maxits = its
                  return
               end if
            else
               eps1=eps*ro
            end if
            eps1=max(eps1,epsmac)
         end if
 
         coeff = 1.0d0 / ro
c        call dscal(n,coeff,ss,1)
         do ii=1,n
           ss(ii,1)=ss(ii,1)*coeff
         enddo
         if (ioutgmr .gt. 0) PRINT *,its,' ',ro
c ** initialize 1-st term  of rhs of Hessenberg system.
         rs(1) = ro
         i = 0
 4       i=i+1
         its = its + 1
         i1 = i + 1
C
C           ss(i1):=A*ss(i)
C
         call ope(n,ss(1,i),ss(1,i1),matvec,va,Ycomm, 
     >            PP, Q, R, Sm, Tm, e1,e5,e2,e4,
     >            UPSTRa,UPSTRb,DOWNSTRa,DOWNSTRb,
     >            rightprecon)

c-----------------------------------------
c  modified gram - schmidt...
c-----------------------------------------
         do j=1, i
            t = ddotp(n, ss(1,j),ss(1,i1))
C           PRINT *,' t0 ',t
            hh(j,i) = t
c           call daxpy(n, -t, ss(1,j), 1, ss(1,i1), 1)
            do ii=1,n
              ss(ii,i1) = ss(ii,i1) - t*ss(ii,j)
            enddo
         end do
         t = ddotp(n, ss(1,i1),ss(1,i1))
C        PRINT *,' t1 ',t
         t = dsqrt(t)
         hh(i1,i) = t
         if (t .eq. 0.d0) goto 58
C        PRINT *,' t ',t
         t = 1.0d0 / t
c        call dscal(n, t, ss(1,i1), 1)
         do ii=1,n
           ss(ii,i1)=ss(ii,i1)*t
         enddo
c--------done with modified gram schmidt and arnoldi step..
c now  update factorization of hh
c---------------------------------------------------------
 58      if (i .eq. 1) goto 121
c-------- perform previous transformations  on i-th column of h
         do k=2,i
            k1 = k-1
            t = hh(k1,i)
            hh(k1,i) = c(k1)*t + s(k1)*hh(k,i)
            hh(k,i) = -s(k1)*t + c(k1)*hh(k,i)
         end do
 121     gam = dsqrt(hh(i,i)**2 + hh(i1,i)**2)
         if (gam .eq. 0.d0) gam = epsmac
C        PRINT *,' gam ',gam
         gam=1.d0/gam
c-----------#  determine next plane rotation  #-------------------
         c(i) = hh(i,i)*gam
         s(i) = hh(i1,i)*gam
         rs(i1) = -s(i)*rs(i)
         rs(i) =  c(i)*rs(i)
c---determine residual norm and test for convergence-
         hh(i,i) = c(i)*hh(i,i) + s(i)*hh(i1,i)
         ro = dabs(rs(i1))
         if (ioutgmr .gt. 0) PRINT *,its,' ',ro
         if ((i. lt. minits) .or. 
     >       (i .lt. im .and. ro .gt. eps1))  goto 4
c
c now compute solution. first solve upper triangular system.
c
C rs := hh(1:i,1:i) ^-1 * rs
 
         call dtrsv('U','N','N',i,hh,maxd1,rs,1)   !+mb
 
c done with back substitution..
c now form linear combination to get solution
         do j=1, i
            t = rs(j)
c           call daxpy(n, t, ss(1,j), 1, sol,1)
            do ii=1,n
              sol(ii) = sol(ii) + t*ss(ii,j)
            enddo
         end do
c restart outer loop  when necessary
      if (ro .gt. eps1 .and. its .lt. maxits) then
        goto 10
      end if
 
 199  format('   its =', i4, ' res. norm =', d20.6)

 
      maxits=its
      eps=eps/eps1*ro !+mb this delivers either actual residual norm or
                      !+mb norm reduction achieved (it depends on iabs)
      return
c------------------------------- end of gmres ----------------------
      end
 
C \end{verbatim} 

      real*8 FUNCTION dnrm2own(N,x,incr)
      implicit none
      integer N,incr,i
      real*8 x(N),sum
C
      sum=0.D0
      DO i=1,n
        sum=sum+x(i)*x(i)
      ENDDO
      dnrm2own=DSQRT(sum)
      RETURN
      END
      
